Examining the evidence for dynamical dark energy 
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We apply a new non-parametric Bayesian method for reconstructing the evolution history of the 
equation-of-state w of dark energy, based on applying a correlated prior for w(z), to a collection of 
cosmological data. We combine the latest supernova (SNLS 3-year or Union2.1), cosmic microwave 
background, redshift space distortion and the baryonic acoustic oscillation measurements (including 
BOSS, WiggleZ and 6dF) and find that the cosmological constant appears consistent with current 
data, but that a dynamical dark energy model which evolves from w < —1 at z ~ 0.25 to to > —1 
at higher redshift is mildly favored. Estimates of the Bayesian evidences show little preference 
between the cosmological constant model and the dynamical model for a range of correlated prior 
choices. Looking towards future data, we find that the best fit models for current data could be well 
distinguished from the ACDM model by observations such as Planck and Euclid-like surveys. 

PACS numbers: 95.36.+X, 98.80.Es 
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Dark energy (DE), the source driving the accelera- 
tion of the universe in the framework of general rela- 
tively (GR), has remained an enigma since its discovery 
[1]. The accumulating observational data, including su- 
pernovae (SN), cosmic microwave background radiation 
(CMB) and large scale structure (LSS), constrains its 
cquation-of-statc w(z), which can be a general function 
of redshift z. An accurate reconstruction of w(z) can help 
us understand DE and gravity: aw/-l would indicate 
a dynamical DE, while a w{z) evolving across —1 would 
imply an additional intrinsic degree of freedom of DE and 
could be a smoking gun of the breakdown of Einstein's 
theory of general relativity on cosmological scales. 

The evolution of w can be reconstructed from data 
using cither parametric or non-parametric methods [2], 
with the latter having the advantage of not assuming any 
ad hoc form of w(z). A common non-parametric approach 
is to bin w in z, or the scale factor a, and fit the bin 
amplitudes to data. This assumes that w(z) is constant 
within each bin, while the neighboring bins are treated as 
independent. But it seems rather unphysical to assume 
perfect correlation of w(z) within a bin, while having 
no correlation between different bins. This approach also 
leads to a practical problem - when fitting many bins to 
data one finds that the uncertainties of binned ui's are 
very large and highly correlated, corresponding to flat 
directions in the likelihood function and a very slow con- 
vergence of Monte Carlo Markov chains (MCMC). Using 
only a few bins improves the convergence, but leads to 
unphysical discrete structures caused by the coarse bin- 
ning. 

A principal component analysis (PCA) shows [3] that 
relatively few uncorrelated linear combinations of the 
bin amplitudes are well-constrained by data, while most 
are practically unconstrained. PCA provides a consistent 



framework for forecasting and comparing the information 
content of future surveys, as well as a tool for data com- 
pression [4]. It has also been used to reconstruct w(z) 
[3, 5] by first performing a Fisher forecast to determine 
the eigenmodes of the covariance matrix based on a par- 
ticular fiducial model, and then fitting the few best con- 
strained modes to data, while setting the amplitudes of 
the poorly constrained modes to zero. This simple trun- 
cation can lead to significant biases and unrcalistically 
small errors in the reconstruction at high z, where the 
data are poor or absent entirely [3]. Ignoring the poorly 
constrained modes is equivalent to the unreasonable as- 
sumption that their amplitudes are measured to be zero 
with infinite accuracy. In addition to biasing, zeroing the 
noisy modes introduces a hidden prior on the smoothness 
of w(z) that is difficult to interpret and quantify. 

To address these issues, we recently introduced a new 
reconstruction method which imposes an explicit prior di- 
rectly on the space of w{z) functions, which is combined 
with observations in a straightforward Bayesian frame- 
work and is simple to implement [4, 6]. (The Gaussian 
Process method [7-9] is close in spirit to our approach, 
but the methods and interpretation of the results dif- 
fer significantly.) The prior is independent of the choice 
of w basis, and can be chosen to prefer smoothness in 
the reconstructions, naturally constraining the flat direc- 
tions in parameter space which are most likely dominated 
by noise. And while some bias in reconstructions is in- 
evitable, particularly when one attempts to reconstruct 
w{z) in regimes where the data are weak, by choosing 
a prior based on theoretically plausible w(z), we ensure 
that the biases are minimised for the most interesting 
models. Here we apply this new method to a collection of 
the latest cosmological data to investigate whether there 
is evidence for dynamical dark energy. 



2 



The choice of the prior is intrinsically subjective be- 
cause it quantifies the initial probability distribution for 
the 'true' model. In principle, it should be determined 
by purely theoretical considerations, e.g. by considering 
ensembles of possible quintessence potentials, but here 
we take a more phenomenological approach. We impose 
a smooth equation of state by assuming w(a) is a Gaus- 
sian field with a covariance described by a correlation 
function 

UI«-^)MHa)-^ fid (a)M«V™ fi V)]) , (i) 

which we take to be translationally invariant in a. We as- 
sume the correlation function to have the form proposed 
by Crittenden, Pogosian and Zhao (CPZ) [4], 

&pz(5a) = Z w (0)/[l + (5a/a c ) 2 ) , (2) 

where a c describes the typical smoothing distance, and 
£ w (0) is a normalisation factor related to a% , the ampli- 
tude of the expected variance of the mean w. The CPZ 
correlation function essentially requires that the w bins 
within the correlation length a c are positively correlated, 
while the correlation falls off when 5a > a c . In [6] we 
have considered other choices of correlation functions and 
found that changing the shape of the correlation function 
does not significantly impact the reconstruction results. 

The two parameters, a c and £ w (0) (or a%) can be 
tuned to adjust the smoothing scale and the strength 
of the prior. A very weak prior (large small a c ) 
results in a noisy reconstruction (large variance), but 
which will match the true model on average (unbiased). 
A stronger prior results in a smaller variance, but pulls 
the reconstructed model towards the peak of the prior, 
which will lead to a bias when the prior conflicts with 
the true model. In [6], we considered reconstructing a 
range of models and found that, by setting a c = 0.06 and 
(7 tB = 0.02, we were able to reconstruct them all without 
a sizeable bias. However, as the functions we considered 
may not fully represent the full range of possible models, 
we also consider a weaker prior with a larger variance in 
the mean, er lB = 0.04, which is still sufficient to curtail 
flat directions in the parameter space. 

We parametrize w(z) using a set of bins w = {wi}, 
i = 1,...,N, spaced uniformly in scale factor a within 
[a m i n ,l]. To make w(z) diffcrcntiablc for the calculation 
of the dark energy perturbations, we interpolate between 
the bins using narrow tanh functions [4]. As long as the 
bin width is small compared to the correlation length, 
the prior largely wipes out the dependence on the choice 
of binning. The number of bins and the range are chosen 
to be large enough so that the results are stable to these 
choices, and alternative binning schemes, e.g. logarithmic 
in z, do not affect the results. Given a binning scheme, 
it is straightforward to construct the prior bin covari- 
ance matrix C [4, 6]. The resulting prior in this space is 

^prior(w) OC e~ (w-w* d ) T Cr 1 (w-w fid )/2 ^ where w fid ig th{J 

fiducial model. 



The reconstructed model is that which maximizes the 
posterior probability, which by Baycs' theorem is pro- 
portional to the likelihood of the data times the prior 
probability, 7>(w|D) oc 7>(D|w) x P prior (w). Effectively, 
the prior results in a new contribution to the total \ 2 
of a model, xprior = (w — w fid ) T C _1 (w — w fid ), which 
penalizes models that are less smooth. One can then use 
standard MCMC techniques to search for the minimum of 
X 2 = Xdata+Xprior- The prior curtails the flat directions in 
the likelihood corresponding to the unconstrained eigen- 
modes of w, giving fast convergence of MCMC chains 
even with a large number of bins. 

We wish to avoid explicit dependence on the fiducial 
model and there are several ways to marginalize over 
it [6]. Here, we adopt the 'floating' average method pro- 
posed in [6], taking the fiducial model to be a flat local 
average of trial bin amplitudes within a range Aa = 0.06. 
Choosing this fiducial model weakens the priors on the 
long wavelength modes, particularly the average mode 
which is given infinite variance. 
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i = 1 


0.06 ±0.2 


0.1 ±0.2 


0.06 ±0.2 


0.2 ±0.2 




i = 2 


-0.3 ±0.3 


-0.3 ±0.3 


-0.2 ±0.3 


-0.2 ±0.4 
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i = 3 


1.3 ±0.6 


0.9 ±0.6 


1.4 ±0.8 


0.8 ±0.8 




i = 4 






0.6 ±0.4 


0.5 ±0.4 



TABLE I. The Ax 2 rows: the improved \ 2 or the best fit 
w(z) model with respect to that of the best fit ACDM model, 
namely, Xw(z) ~ Xacdm- The S/N row: the significance of 
w(z) deviating from —1, which is simply \/|Axtotl- The en 
rows: the mean and 68% CL error on the principal component 
amplitudes. 

We apply our method to a joint dataset of the latest 
cosmological observations including SN, CMB, LSS and 
the measurement of Hubble parameter at various red- 
shifts. We compare two different SN samples: SNLS 3- 
year [10] and Union2.1 [11]. We used the WMAP 7-year 
CMB spectra [12], the H(z) data compiled by [13], the 
/<78 estimates from the redshift space distortions (RSD) 
and the baryonic acoustic oscillations (BAO) measure- 
ments from SDSS-II [14], SDSS-III BOSS [15], 6dF [16] 
and the WiggleZ survey [17]. Given this joint dataset, 
we use MCMC [18] to sample the parameter space P = 
(u>b,u} c ,@ s ,T,n s ,A s ,w-L,...,W2o,J^) where uj b and w c are 
the baryon and cold dark matter densities, O s is the ra- 
tio of the sound horizon to the angular diameter dis- 
tance at decoupling, r is the optical depth, n s and A s 
are the primordial power spectral index and amplitude 
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FIG. 1. Panels (A1-A4): The best fit w(z) with marginalized 68 and 95% CL error (shaded bands) reconstructed from a joint 
data set of the latest data including SN, CMB, H(z), RSD and BAO. This best-fit w(z) is used to construct mock data for 
future Planck and Euclid-like surveys, and the resulting reconstruction is shown by the magenta dots (with 68% CL error 
bars.) The four panels show the result using different combinations of priors (standard and weak) and supernovae data sets 
(SNLS3 and Union2.1). Panels (Bf-B4): The 68% CL constraints on w(z) using different data subsets: SN+CMB+#(z)(black 
solid), SN+CMB+#(2)+RSD(magenta dashed), SN+CMB+i/(z)+RSD+BAO(blue band). Panels (C1,C2): These recon- 
structed models affect the fit to the BAO and RSD data; we compare the reconstruction from the combined data (using 
SNLS3 and weak prior) with the ACDM model (green dashed). 



and w\, ...,u>2o denote the 20 w bins uniform in a from 
flmin = 0.4 to a = 1, corresponding to z G [0, 1.5]. We fix 
w = — 1 at z > 1 .5 1 . We also marginalize over parameters 
accounting for the calibration uncertainty in measuring 
the intrinsic SN luminosity. We use a modified version of 
CAMB [19] to calculate the observables and include dark 
energy perturbations for an arbitrary w following the pre- 
scription in [20]. 

Fig. 1 shows the reconstruction results obtained with 
MCMC sampling of the parameter space P to find the 
minimum of x^ata + Xprior- The shaded regions in panels 
(A1-A4) illustrate the 68 and 95% CL uncertainty of the 
reconstruction, while the solid curves show the best fit 
w(z) models. We show separately the results using stan- 
dard or weak priors, and using the SN data of SNLS3 
or Union2.1, combined with all the other data. Looking 
at the different priors, we see similar trends, but the re- 
construction uncertainties for the weaker prior are larger 
at low z(z ~ 0) and high z(z > 0.6), where the con- 
straints from data are weak. At the sweet spot, z ~ 0.15, 
the error (A«;| z= o.is — 0.12) is independent of the prior, 
since there the data are sufficiently strong to outweigh 
the priors. Interestingly, the most likely reconstruction 
using the SNLS3 data favors a transition from w < — 1 



Adding a w bin from z = 1.5 to 1100 has a negligible effect on 
our results. 



at low redshift to w > — 1 at higher redshift, a behav- 
ior that is consistent with the quintom model which al- 
lows w to cross —1 [21]. The result for Union2.1 shows 
a similar trend, but at a lower significance. As shown in 
Table I, for SNLS3 combined with all other data, the sig- 
nificance of w departing from —1 is 2.3(2.5)-cr under the 
standard(weak) prior. 

To understand this result, we show the reconstruc- 
tion ranges using different data combinations in panels 
(B1-B4). Without the RSD and BAO data (black solid), 
SNLS3 data slightly prefers a model with w < — 1 at 
z ~ 0.25, but Union2.1 is very consistent with w = — 1. 
With the RSD data (magenta dashed), the reconstruction 
at z > 0.25 is better constrained, and w(z) is slightly in- 
creased at z > 0.25. After adding BAO, the evolution 
from w < —1 at low redshifts to w > — 1 at high red- 
shifts becomes more evident, particularly when using the 
SNLS3 data; this suggests that the BAO data is at least 
partially responsible for the trend. As shown in panel CI 
for the SNLS3 data and a weak prior, the best fit w(z) 
can fit the BAO data (the comoving angular diameter 
distance Dy) better than the ACDM model, especially 
for the BOSS BAO data point at z ~ 0.57 (the x 2 for 
BAO is reduced by 1.8 as shown in Table I). This natu- 
rally explains the result we found: w needs to be less than 
— 1 at low redshifts to get a higher Dy at z ~ 0.6 com- 
pared to ACDM, while CMB requires the average w to 
be close to —1, so the high redshift behaviour is pushed 
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FIG. 2. Panels (A1,A2): Eigenvalues of the covariance matrix of the w bins obtained from the data plus prior (diamonds 
and circles) and prior alone (shaded region) using MCMC. Panels (B1,B2): The first 3 (4) data-dominated eigenmodes of the 
covariance matrices. These panels show only small differences arise when different SN datasets (SNLS3 and Union2.1) are 
combined with other data. Results for the standard prior (Af ,B1) and weak prior (A2,B2) are shown separately. Panels CI 
and C2 show the improved likelihood (AlnL), the reduced fraction of the sampled parameter space (AlnV), and the logarithm 
of the evidence ratio as a function of <7bm relative to the ACDM model (which is the small ahm limit.) Cf and C2 are for the 
present data (SNLS3) and the simulated Planck and Euclid-like data respectively. 



above —1 to compensate. This quintom behavior is also 
favored by SNLS3, where the contribution to \ 2 from SN 
is reduced by 2.2. 

We also perform a PCA on the posterior distribution, 
identifying the uncorrelated collective w(z) degrees of 
freedom after marginalizing over the other parameters. 
In Fig, 2, panels (A1,A2) show the eigenvalues of data 
plus prior (points) and prior alone (shaded regions) . The 
choice of the standard (weak) prior allows the informa- 
tion in the strongest 2(3) data modes to survive the 
addition of the prior. In panels (B1,B2), we show the 
data-dominated eigenmodes ei(z). We can project w(z) 
onto this orthonormal basis, i.e. 1 + w(z) = J^. oiiei{z). 
Note that if w = — 1, all the a's should be zero. Pro- 
jecting the best fit w(z) models onto the eigenbasis, 
we obtain the constraint (central value and 68% CL 
error) on the a's shown in Table I. We can see that 
A% 2 ~ X)i=i[ Q; i/ <T ( a i)] 2 j meaning that the improvement 
in the fits is almost entirely accounted for by these three 
or four principal components. 

The reconstructed models necessarily provide a better 
fit than ACDM (Table I) , but the crucial issue is whether 
the improvement in fit is enough to compensate for the 
increased parameter space in the prior model. Naively, 
the non-parametric binning has 20 additional degrees of 
freedom, but this far overestimates the true parameter 
volume when the correlated prior is imposed; in the limit 
of a very strong prior, the model can become equivalent 
to the cosmological constant model itself. 



To quantify this, we estimate the evidence ratio, called 
the Bayes' factor, within a family of models which inter- 
polates smoothly between the prior described above and 
ACDM; this is implemented by adding a larger and larger 
diagonal term to the prior which shrinks the variance in 
each bin, X prior = (w + 1) T [C- 1 + a^I] (w + 1), and 
shifts all the eigenvalues by a constant. We calculate the 
evidence [22], E = J d"P P(D|P) P pl . ior (P), by assum- 
ing a Gaussian prior for all parameters P, described by a 
more general covariance matrix, C pr i or , and also assume 
the resulting posterior is a Gaussian function of the pa- 
rameters around the best fit, with covariance described 
by C post . We find: 

Eoc( detCpoat ) 1/2 c~^' /2 , (3) 

V det Cp r ior / 

where Xh f ^ s the minimum total x 2 given the data and 
prior. The first term quantifies the fraction of the initial 
prior parameter volume which remains consistent with 
the measured data, while the second term accounts for 
how well the model is capable of fitting the data. 

We plot the logarithm of the Bayes' factor in Fig. 
2, comparing a range of models approaching ACDM 
to ACDM itself. For the original prior discussed above, 
marginalizing over fiducial models formally set the vari- 
ance in the mean to infinity, meaning detC pr ior diverges 
and the evidence ratio approaches zero. However, choos- 
ing a more reasonable variance in the mean equation of 
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state, one finds evidences comparable to ACDM. Perhaps 
surpisingly, for some choices of this variance, dynamical 
dark energy is actually preferred over ACDM. However, 
choosing a particular prior after seeing the evidence is 
post hoc, and to be convinced one would need to see it 
for a range of priors most people would find reasonable. 

While current data is unable to distinguish ACDM 
from the best fit models, it is interesting to see what the 
future data could do. For this purpose, we perform a fu- 
ture forecast assuming that today's best fit model was the 
true model. We assume a future dataset of a deep SN sur- 
vey [23] which is similar to the expected Euclid mission 
[24], and the CMB data from Planck[25]. The reconstruc- 
tion can be performed using the Wiener filter projection 
[6] and the reconstruction is shown in Fig 1. The yellow 
dots with error bars show w(z), which would be recon- 
structed with very little bias. In both cases, the dynam- 
ical models can be clearly distinguished from w = — 1. 
This is supported by the evidence calculation in the lower 
right of Fig. 2; for the Euclid-like data, the evidence ratio 
is large and sustained over a wide range of priors. 

Our non-parametric reconstruction is simple to imple- 
ment and allows for straightforward interpretation and 
calculation of evidence ratios. While we see tantalizing 
hints for dynamical dark energy, the case is by no means 
conclusive. It is interesting that, except for Union2.1, the 
fits to all the data are improved with our reconstructed 
model (Table I). Note that while SNLS3 and Union2.1 
overlap significantly, SNLS3 has a more homogeneous 
higher redshift sample, which may make it less suscepti- 
ble to systematic errors. Future data will be able to dis- 
tinguish this reconstruction from ACDM, while Planck 
and the full three-year SDSS SN sample might help in 
addressing this question in the near term. 
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